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i ABSTRACT 

Numerical  methods  are  presented  for  determining  the  stability  of  a 
multiple  autoregressive  scheme.  The  question  of  stability  is  important  in 
prediction  theory,  but,  furthermore,  the  methods  for  determining  stability 
are  such  that  they  yield  information  which  is  important  when  fitting  time  series 


models  to  data. 
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1,  Introduction 

Let  [X(t),teZ)  , Z the  sot  of  a.11  integers,  be  a zero  mean, 

d-dimensional,  covariance  stationary  time  series  with  absolutely  summablo 

T T 

autocovariance  function  R(v)  c E(X(t)X  (t  + v))  , vcZ  , A denotes  the 
transpose  of  the  matrix  A . Then  the  spectral  density  matrix, 

1 * 

f(i«)  = 2;;  S'  R(v)  , oi)e[-w,ir] 

V s . «p 

and  the  R(«)  , are  Fourier  pairs, 

f(UJ)  e^'^*  dtt) , veZ  . 

Formally  define  the  inverse  spectral  density  matrix  by 

1 -♦ 

fi(uu)  s — 2^  • UDe[-w,  ^]  » 

4v 

e 

where  A denotes  the  complex  conjugate  transpose  of  the  matrix  A . If  tlisre 
exist  positive  finite  constants  X2  such  that 

Xjl^  < £(uu)  < X^I^  . (1.1) 

where  1.  denotes  the  d-dimensional  identity  matrix,  then  fi(u))  exists  and 

G 

has  integrable  components  and  we  can  define  the  covinverse  matrices 


Ri(v)  ■ f*  fi(ttJ)  dw  , veZ 

• •Tf 
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Further,  (Masani  (1966))  there  exirits  an  infinite  order  autoregressive  repre- 
sentation for  {X(t),  teZ  J 


Z A (j)X(t  - j)  = e(t)  . (1.  2) 

^ CD  ^ 

j=0 


where 


e(t)  = X(t)  - E(X(t)lX(t  - 1)  , X(t  - 2), . . . ),  teZ  , 

is  the  infinite  memory  one-step  ahead  prediction  error  and  is  a multiple  white 
noise  series,  with  positive  definite  covariance  matrix  Z , 

E(e(t)/(t  t v)»  . ^ , 

where  6.  ,=l,ifj=k,  and  zero  otherwise.  Defining  L to  be  the  back- 
Jt  * 

ward  lag  operator,  L^X  (t)  = X(t  - j)  , jeZ  , we  can  write  (1.  2)  as 


G (L)X(t)  = c(t)  , 

CO  ^ ^ 

where  the  complex  valued  matrix  polynomial  G^(«)  is  defined  by 

CO 

G (z)  = Z A 0)*^  . 

w CO 
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The  representation  (1.  2)  affords  a representation  of  f(*)  : 


f (U))  r G"  V“^  , tt)€{-ir,ir]  . 

00  iT  00  00  00 


Definition  of  Autoregressive  Approximant 
Let 


fp(uj)  = ^ G‘\e‘  ‘^rpG‘*{e*‘^  , tt)e(-iT,  ir] 


where 


G„(«)  . Z A O)*'’  . 

P j.O  P 


and  A (1),...  ,A  (p),  2 are  solutions  to  the  Yule-Walker  equations, 
P P p 


Z A (j)R(j  - v)  = 6 Z , v=0,...,p. 

j=0  ^ v,u  p 


(1.3) 


We  call  f (•)  the  p^  order  autoregressive  approximant  to  f(*)  (Parzen 
P 

(1975)).  The  use  of  the  term  autoregressive  is  natural  since  the  p-memory 
linear  prediction  error,  , defined  by 


(t)  = X (t)  - E(X  (t)  |X(t  - 1), . . . , X(t  - p)) 

M ^ ^ ^ 


is  given  by 
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= S A^(j)X(t  - j)  , A (0)  = I , (1. 4) 

j=o  ^ 

where  A (1) A (p)  . and  Z = E(e^P^t)c^P^'^(t)) 

P P p /w 

are  obtained  from  (1.3). 

Looking  at  autoregressive  approximants  from  a prediction  theory  point 
of  view  we  see  that  for  stability  of  predictions,  it  is  desirable  that  we  be  ?bJe 
to  write  X (t)  as  a linear  combination  only  of  past  and  present  : 

09 

X(t)  = Z B (k)e^P^t  - k)  , teZ  . (1. 5) 

^ D ^ 09 

k=0  P* 

Definition  of  Stationary  Autoregressive  Approximant 

An  autoregressive  approximant  of  the  form  (1.4)  which  has  a represen- 
tation, as  a limit  in  the  mean,  of  the  form  (1.  5)  is  said  to  be  stationary. 

The  use  of  autoregressive  approximants  for  d = 1 has  become  quite 
widespread  (Parzen  (1974),  Akaike  (1971)).  Parzen  (1976)  has  recently 
extended  his  method  of  determining  the  order  of  autoregressive  approximants  to 
the  multiple  case.  In  Section  2 we  develop  theory  and  algorithms  for 
determining  whether  a given  set  of  A(l), . . . , A(p)  , and  Z determine  a 
stationary  autoregressive  approximant.  This  theory  allows  us  to  obtain  two 
further  results: 

(a)  Autoregressive  approximants  for  all  orders  from  a covariance 
stationary  time  series  whose  spectral  density  is  bounded  f bcve 
and  below  by  positive  finite  constants  are  stationary. 


6 


(b)  Algorithms  for  determining  the  autocovariances  R(0), . . . , R(p) 
corresponding  to  the  A(*)  and  £ of  a stationary  autoregres- 
sive approximant.  These  algorithms  have  several  applications 
in  multiple  time  series  analysis,  two  of  which  are  given  in 
Section  3. 

Finally,  in  Section  4 we  give  a numerical  example  of  the  algoritlims 
developed  in  Section  2. 

2.  The  Algorithms 

Given  matrices  A(l), . , . , A(p)  and  £,  define  the  autoregressive 
process  [X(t),teZ)  by 

P 

S A(j)X(t  - j)  * c(t)  , (2.1) 

j=0 

where  {je(t),  teZ)  is  a white  noise  process  having  covariance  matrix  S , 

We  develope  algorithms  for  determining  whether  there  exists  a stationary 
solution  to  (2.  1)  in  terms  of  only  present  and  past  e(*)  . 

Lemma  1;  (Hannan  (1970,  p.  14)) 

If  (£(t),  teZ}  is  a white  noise  series,  then  given  p initial  values, 
X(-p), . . . , X(-l)  , the  solution  to  (2.  1)  is  uniquely  defined.  If  all  zeros  of  the 
determinant  | G(z)  | are  outside  the  unit  circle,  the  solution  is  of  the  form 

00 

X(t)  = Z B(k)cft  - k)  , teZ  , 
k=0  ~ 

where  the  elements  B^j(k)  of  the  B(lc)  converge  to  zero  exponentially  as 


k -•  * . 
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The  B(«)  can  be  obtained  by  comparing  coefficients  of  z*'  in 


G’\z)  = H(z)=  r B(k)z^  . 

k=0 


Method  1 (via  the  Schur  matrix) 

Method  1 consists  of  determining  whether  the  zeros  of  |g(z)|  are 
outside  the  unit  circle.  Writing  |g(z)|  as  a polynomial  of  degree  pd  , 


pd 

|g(z)|  = Z . 

j=0 


(see  Robinson  (1967),  p,  162  for  example)  we  can  implement 


Theorem  on  the  Schur  Matrix  (Schur  (1917)) 

n 

The  zeros  of  the  polynomial  g(z)  s 2 a(j)z^  are  outside  the  unit 

j=0 

circle  if  and  only  if  the  matrix  S defined  by  S = (S.,  ) , 


V 

S.,  = S (a(j  - 1 - 4)a(k-  1 - 4)  - a(n+  1 +■(.- j)a(n+  1 +‘t-k)J  , 
4=0 


jfk  = l,...,n,  v = min  (j  - l,k  - 1)  , is  positive  definite.  The  matrix  S ic 
called  the  Schur  matrix  of  the  polynomial  g(«)  . It  is  particularly  useful 
when  d = 1 , since  then  (Pagano  (1973)),  if  the  zeros  of  |G(z)j  are  all  outsid 
the  unit  circle, 


2 

a 


r 

p 


9 


(^?) 
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where  is  the  p x p Toeplity  matrix  having  R(j  - k)  in  the  row  and 
, th  , . 2 2 

K column,  while  a = E(e  (t))  , Unfortunately  there  appears  to  be  no  vector 
analogue  of  (2.  2). 

To  decide  whether  S is  positive  definite  one  may  calculate  its  modified 
Cholesky  decomposition  (see  Wilkinson  (1967)).  Then  S is  positive  defir.:fa 
if,  .and  only.  If,  the  diagonal  elements  of  the  diagonal  matrix  in  its  modified 
Cholesky  decomposition  are  positive. 

As  mentioned  above,  method  1 for  d = 1 produces  the  autocovariancf'.T 
* '^'bile  for  d > 1 it  does  not.  The  methods  given  below  do  give  th? 
R(*)  , They  are  based  on  the  following  lemmas. 

Lemma  2 (Hannan  (1970),  p.  19) 

Let  R(v)  = E(X(t)X^(t  + v))  , veZ  . If  the  zeros  of  |g(z){  are  all 
outside  the  unit  circle,  then  the  R(*)  se<iuence  satisfies 

P 

A(j)R(j  - v)  s 6 L , V a 0 . (2.  3) 

j=0  ° 


Lemma  3 

If  Z is  positive  definite,  then  the  zeros  of  |g(z)|  are  all  outside 

the  unit  circle  if,  and  only  if,  F is  positive  definite,  where  F is  the 

P p 

p X p block  Toeplitz  matrix  having  R(j  . k)  in  the  row  and  k*^^  colurv'. 
of  blocks,  i.  e.  , F^  = Toepi(R(0), . . . , R(1  - p)  . 
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Proof 

Hannan  (1970),  p,  329,  shows  sufficiency,  while  Whittle  (1963)  shows 
necessity. 

Theorem;  Bounded  Spectra  Guarantees  Statiox^ry  Approxlmacts 

If  [X(t),  teZ)  is  a covariance  stationary  time  series  whose  spectral 
density  matrix  satisfies  (1.  1),  then  the  autoregressive  approximants  for 
p ^ 1 defined  by  (1.4)  are  stationary. 

Proof;  . 

We  show  that  Z and  F are  positive  definite  for  p ^ 1 , thus 
^ P ^ 

ensuring  that  the  zeros  of  | Z A (j)z^  | are  outside  the  unit  circle. 

j=0  P 

X X • X 

Define  the  nonnull,  pd-vector  K>  =('t,I.,.I't'  ),  where 

~ ~1  . .~p 

T 

. is  a d-vector.  Then  T is  positive  definite  if,  and  only  if,  6 F v > 0 

p ' ' ~ p/s. 

for  all  such  K . Now 

T P T 

I T 1=  Z 6.  R(j  - k)6, 

~ P~  j,k=l~J 


P 

= Z 


j»k=l 


^i(j-k)U) 


* C*(u))f(u})C((»)  du? , 


where 


P 

C(uo)  = Z 

j=l 


is  nonnull.  Now  (1.  1)  implies 
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-t^r  t s X-  c*(uj)c(uj)  d 

~ p~  1 J_ir  ~ ' 


uu 


= X. 


uu 


= X 


1 


T 

■t  . > 0 . 


Thus  we  can  write  the  Yule-Walker  equations  (1.  3)  as 


FA  = -R  , 
P P P 


Z = R(0)  - R^r"^R  , 
P P P P 


where  A = (A  (1)  I ^ A (p))*^  , and  R = (r"^(1)  I ...  I R (p))  . 

p p . . p P . . 

r and  r are  positive  definite,  and  by  the  clockwise  rule  of  matrix 

P P+1 

determinants. 


R(0) 

R 

P 


r 

p 


= |r  I |r(0)  - R^r'V  I 

p p p p 


= r 


1^1 


Since  Z is  a covariance  matrix,  the  proof  is  complete. 
P 


Nov/ 
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Method  2 (via  the  Yule- Walker  equations) 

Method  2 consists  of  solving  (2.  3)  for  the  R(»)  and  then  checking  the 

resulting  for  positive  definiteness.  The  transformation  from  the  A(®)  to 

the  R(«)  is  not  obvious  (see  Quenoville  (1957)).  The  symmetry  of  R(0)  and 

the  noncommutativity  of  the  R(»)  and  A(»)  seems  to  rule  out  a solution  of 

(2,  3)  in  terms  of  the  matrices  themselves.  The  matrix  notation  in  (2.  3) 

2 

disguises  the  fact  that  there  are  in  fact  (p  + l)d  scalar  equations  in  the  fire- 

p + 1 matrix  equations.  Further,  the  symmetry  of  R(0)  and  E actually 
2 

leaves  s = pd  + d(d  + l)/2  distinct  scalar  knowns  and  s unknowns  in  these 
equations,  namely 


P 

2 


d 

E 


j=0  tsl 


A ,(j)R.  (j  - V)  = 6 -Z 
m-l'  tn  V,  0 mn 


(2.  41 


for  V s 0,«a.,p  s m,  n = l,...,d. 

Define  the  vec  operator,  vec  (C)  , on  the  (n  X m)  matrix 

C = (C  ,..,,C  ) by 

~ 1 ~ m 


vec  (C)  = (cf  : ...  : )^ 


and  the  Ivec  operator,  Ivec  (C)  , as  the  vec  operator  on  the  lower  triangulai 
portion  of  C . Then  define  the  (s  x 1)  vectors  r and  g by 


r = (r^(0)  : . ..  : r^(p))^ 


& = (&^(o):  ... : g'^(p))'^ , 


• • 
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where 


J Ivec  (R(0))  , j = 0 

' vec  (R(j»  f j = 1. . • • tP 


&<j)  = 


Ivec  (S) 


j = 0 


where  0 denotes  the  n-vector  of  zeros.  Then  the  equations  (2.  4)  can  be 
~n 

written 


Wr 


= & 


* 


(2.5) 


where  W is  an  (s  X s)  matrix  determined  as  follows. 

T 

Since  R(-v)  = R (v)  , we  have 

R^j^(v)  , v=l 

R.^(v)  = / , v=-p,...,-l  (Z.':) 

R (0)  , V = 0 , 

m m_ 

^ 1 2 


where  m^  = max 
is  given  by 


m^  = min  (j,k) 


The  row  number  of 


R.k(v)  in 


i . (V)  = 


(k  - Id  + j - k(k  - l)/2  , V = 0 

d(d  + l)/2  + (v  - l)d^  + (k-l)d+j  , v=l,...,p 

j»k  — l|...|d 
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Thus  to  determine  W : 

(a)  Set  W initially  equal  to  zero 

(b)  For  fixed  m , n , v , j , t in  (2. 4)  : 

(i)  The  row  number  of  W corresponding  to  m , n , v 

is  given  by  t,  = i (v)  . 

1 mn 

(ii)  Find  by  (2.  6)  the  indices  V , n’  , u'  of  the  represen- 
tation R.,  (u')  of  R,  (u)  , u = j - V , for  nonnegative 

u'  . 

(iii)  Add  A (j)  to  W where  t_  = i.  (u')  . 

m'C'-’  tjt^  2 t'n' 

If  (2,  5)  cannot  be  solved  uniquely  for  r , then  the  process  is  not 

stationary.  If  (2,  5)  can  be  solved,  it  is  not  difficult  to  reorder  the  elements 

of  r to  form  R(p), , . . ,R(p)  , and  T , and  to  check  P for  positive 
~ P P 

definiteness. 

Note  that  it  is  possible  to  solve  (2,  5)  uniquely  for  r while  P is  not 

~ P 

positive  definite.  For  example,  let  p = d = 1 , Then  (2.  5)  becomes 


If  |a|  > 1 and  > 0 , then  P^  = R(0)  = a^/(l  - aS  < 0 . 
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Method  3 (iterative) 

If  d is  large,  the  value  of  s may  prohibit  the  use  of  Method  2 unless 
one  uses  an  equation  solver  to  solve  (2,  5)  which  does  not  require  storage  of 
the  entire  W matrix.  Thus  we  propose  an  iterative  method  for  calculating 
which  requires  much  leas  storage  than  the  Yule- Walker  equation  method. 

By  lemma  1,  if  the  zeros  of  |g(z)|  are  all  outside  the  unit  circle, 
we  can  write 


G'^z)  = H(z)  = Z B(k)z^  , (2.7) 

ksO 

where  the  (“t, m)  element  of  B(k)  converges  exponentially  to  zero  as 
k - ® , i.  e. , there  exist  constants  C 0 , p such  that  for  k sufficiently 
large. 


• |c|  <».  IpI  <1. 


Thus  we  can  obtain  R(0), . , . ,R(p)  by 


R(v)=  Z B(k)ZBV  + v)  , (2.8) 

k=0 

The  B(»)  can  be  obtained  by  solving  the  block  triangular  system  of  equations 


G(z)H(z)  . . 


(2.  9) 
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Equating  coefficients  of  in  (2.  9),  we  have  the  difference  equation  with 
B(0)  = . 


minf-t,  p) 

r A(j)B(t  - j)  = 0 , l>  0 . (2.  10) 

j=0 


R(v)  can  be  approximated  by 


R (v)  = S B(k)2B  (k+v)  , v = 0 p.  (2.11) 

^ k=0 


Then  we  have  the  following. 


Convergence  Theorem  for  R.  .(•) 

M 


Let  lie  11  denote  a norm  of  an  (n  X n)  matrix  C = (C  ),  Then  a 
necessary  condition  for  A(l), . . . , A(p)  , and  Z to  define  a stationary  auto- 
regressive approximant  of  order  p is  that 


llR(v)  - Rj^(v)ll  = 0(y'^>  , 1y|  < 1 , 


uniformly  In  v = 0,  ,..,p. 

Proof; 

From  (2.  8)  and  (2. 11)  we  have,  if  the  zeros  of  |g(z)|  are  all  outside 
the  unit  circle. 


I 


Thus  the  iterative  method  consists  of  the  following  steps: 

(a)  Find  B(0), . . . , B(p)  and  Rp(0), . . . ,Rp(p)  by  (2.  10)  and  (2.  1 1 ). 

(b)  Given  B(M  - p), . . . , B(M)  and  R^JO) R^,(p)  , find 

M M 

B(M  + 1)  by  (2, 10)  and  the  i’Y 

(c)  Iteration  stops  if  for  a given  6 , 


17 


If  the  algorithm  does  not  converge,  the  approximant  is  not  stationary. 

If  convergence  is  achieved,  can  be  checked  for  positive  definiteness. 

Comparison  of  the  Methods 

(a)  In  most  cases  knowledge  of  the  R(*)  is  desirable,  so  methods 
2 and  3 seem  preferable. 

(b)  To  write  a computer  program  of  method  2 requires  approximate!; 

2 4 2 

p d storage  locations,  while  method  3 requires  3(p  + l)d 

locations  (A(l ),...,  A(p),  Z,  B(M  - p), ...,  B(M),  and 

Rj^(O), . . . • 

(c)  Preliminary  empirical  evidence  indicates  that  the  value  of  M 
necessary  to  attain  convergence  in  method  3 is  not  large. 

(d)  Method  3 does  not  have  the  problem  of  the  matrix  W of 
method  2 being  ill-conditioned.  However,  method  3 is  an 
iterative  method  and  may  suffer  the  difficulties  of  any  iterative 
algorithm. 

Thus,  unless  the  values  of  p and  d are  quite  small,  method  3 appears 
to  be  the  most  useful  algorithm. 


Applications  of  the  Algorithms  to  Obtain  R(*) 


An  Application  to  Scalar  Mixed  Time  Series 

Consider  the  scalar  mixed  autoregressive  moving  average  process, 
[X(t),teZ}  , of  order  (p,  q)  , defined  by 
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P *1 

L a(j)X(t  - j)  = r 0(k)e(t  - k)  , teZ  , 
j=0  k=0 


where  {e(t),teZ}  is  a series  of  uncorrelated  random  variables  with  zero 

2 

mean  and  variance  a , and  the  zeros  of  the  complex  polynomials 


g(z)  = Z a(j)z'^ 
j=0 


h(a)  = r 3(k)z 
k=0 


are  outside  the  unit  circle.  The  information  matrix  (Whittle  (1953))  of  the 

T T 

parameters  a = (a(l ),...,  a(p))  and  8 = (3{1)» , . . , 3(9))  is  given  by  (see 
Hannan  (1970),  p.  392) 


1 = 


,aa 


o3.T 


d'^) 


,33 


where 


, i(j  - k)u} 


iU).  i2 


j#k*l,.,.jp  • 
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i(j-k)o) 


(r''),  =4-  r -= — ^ — T dll)  , j,k  = i,...,q  , 


^ 1 i(i-k)u) 

aS  1 (.w  e , 

(I  1 : 1 — dyj  . J = 1»....P 

' 'ik  2ir  J.w  iu;  -iu), 

J g(e  )h(e  ) 

k — i9»«c>q 


Consider  the  two  dimensional  autoregressive  process  u(»)  of  order 
m = max  (p,q)  defined  by 


where 


G(L)u(t)  = e(t)  , 


u<t)  = (u,(t)  , a,{t))^  , e(t)  » (e(t),  , 

^ i c ^ 


m ra'O)  0 


G(a)  = Z 

j=0  0 3’(j) 


g(*)  0 
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where  a'(j)  = a(j)  for  j = 0,...,p  , and  zero  otherwise,  and  3'(j)  = P(j) 
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= Ryjl  jO  - It)  . = 1,.  . . ,P  , (3.1) 

33 

^jk  = ^‘u;2,2<^  ■ ’ j,k  = l,....q,  (3.2) 

(I^)jk  = ' j = l,...,p.  k = l,...,q  . (3.3) 


Thus  to  find  I in  a convenient  way,  without  any  numerical  quadrature,  we 
can  find  the  R^(*)  hy  the  methods  of  Section  2,  and  form  I from  (3.  1),  (3.  2), 
(3.  3).  This  is  of  particular  importance  when  estimating  the  parameters  of 
the  scalar  mixed  process  (see  Hannan  (1970),  p.  390). 


An  Application  to  Constrained  Autoregressive  Estimation 


Given  a sample  x(l ),...,  x(T)  , consider  the  problem  of  efficiently 
estimating  the  autocovariance  sequence  R(*)  of  a multiple  autoregressive 

process  X(»)  where  the  coefficients  a = vec  (A(l)  I ...  lA(p))  are  subject 

* • • 

to  the  linear  constraints 


Ba,  = a,  . 

^4  1 


(3.  4) 


Hannan  and  Terrell  (1972)  consider  problems  of  a similar  nature. 

A A A 

Efficient  estimates  A(l), . . . , A(p)  , and  Z of  A(l),. . . , A(p)  and  Z 
(ignoring  the  constraints)  are  obtained  by  solving  the  Yule-Walker  equations 
(2.  3)  with 
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where  T = TOEPL{R_(0).  . . . , R_(l  - p)  3 and 

P ^ ^ 

C ® D denotes  the  Kronecker  product  of  the  (n  X m)  matrix  C = 
the  (r  X s)  matrix  D = 

C ®D  = 

Then  an  efficient  estimator  a,  of  a,  is  given,  by  applying  the  theory  of 
optimisation  under  linear  constraints  (Pagano  (1974)),  as  the  solution  of  the 
normal  equations 


T 

(B^  I 


Then  to  obtain  efficient  estimators  R(l)  ,...,R(p)  of  R(0), R(p)  subject 
to  (3.  4),  one  can  use  the  methods  of  Section  E, 
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4.  Numerical  Example 

Let  p = 2,  d 5=  2 , and 


-.  2580 

1429 

» 

A(2)  = 

1471 

. 1280 

. 5572 

-. 8637 

» 

4732 

. 3305 

2.  6034  .9053 


. 9053  2.  1450 


Method  2 gives  the  exact  values  of  R(») 


After  9 iterations,  method  3 gives  values,  correct  to  6 = 1 x 10 


R(0)  = 


2.9999974  1.0000470 

1.0000470  4.0000683 
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R(l)  = 


I 


i .9999955 
. 5000751 


-.  4999354 
2.  0000473 


R(2)  = 


.499851i 

.0499245 


. 1000909 
. 5999987 


I 


Thus  for  this  second  order  autoregressive  process,  one  needs  an  approxi- 
mating moving  average  process  of  order  11  to  obtain  five  place  accuracy. 
(Calculations  performed  on  a CDC-6400.  ) 


I 


f 


I 
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